Hydrodynamic Correlations slow down Crystallization of Soft Colloids 
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Crystallization is often assumed to be a quasi-static process that is unaffected by details of particle 
transport other than the bulk diffusion coefficient. Therefore colloidal suspensions are frequently 
argued to be an ideal toy model for experimentally more difficult systems such as metal melts. In this 
letter, we want to challenge this assumption. To this aim, we have considered molecular dynamics 
simulations of the crystallization in a suspension of Yukawa-type colloids. In order to investigate 
the role of hydrodynamic interactions (His) mediated by the solvent, we modeled the solvent both 
implicitly and explicitly, using Langevin dynamics and the fluctuating Lattice Boltzmann method, 
respectively. Our simulations show a dramatic reduction of the crystal growth velocity due to His 
even at moderate hydrodynamic coupling. A detailed analysis shows that this slowdown is due 
to the wall-like properties of the crystal surface, which reduces the colloidal diffusion towards the 
crystal surface by hydrodynamic screening. Crystallization in suspensions therefore differs strongly 
from pure melts, making them less useful as a toy model than previously thought. 
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Crystallization and nucleation of undercooled melts 
are often studied using model systems of charged col- 
loids in solution [lfj, |l7J, |l8|, such as polystyrene (PS) 
or polymethylmethacrylat (PMMA) spheres suspended 
in water [7J, [la, |26j . The solvent gives rise to hydrody- 
namic interactions (His) between the colloidal particles 
that are not present in, e.g., metal melts. The influence 
of His on the dynamical properties of colloidal suspen- 
sions has been extensively studied in recent years [231.124 1. 
Lowen et al. [171, |l9( showed that the ratio of the long- 
time to short-time self-diffusion coefficients has a univer- 
sal value along the fluid freezing line. Recent studies by 
Pesche [25[ and Nagele [20j of quasi-2D dispersions show 
that His have an impact on the self-diffusion function 
in these soft-sphere suspensions. However, since crystal 
growth happens on much larger time scales, it is com- 
monly believed to be a quasi-static process that is unaf- 
fected by His. For example, in classical nucleation theory 
(CNT) [14j, hydrodynamic interactions are assumed to 
enter only through the effective diffusion constant of the 
attaching colloids, which can be measured conveniently 
in the bulk liquid. Also in most computer simulation 
studies of nucleation, hydrodynamic interactions are ne- 
glected to avoid the high computational costs. 

In the following, we will show with the help of com- 
puter simulations that His do have a remarkable influence 
on the dynamics of the crystallization in a colloidal sus- 
pension. Even at moderate coupling, we found the crystal 
growth velocity to be reduced by a factor of three. Sim- 
ilar findings have been reported by Schilling et al. using 
different simulations techniques 27] . 

To investigate the influence of hydrodynamic interac- 
tions on crystal growth, we studied the crystallization 
of Yukawa-type particles confined between two planar 
walls. The influence of hydrodynamic correlations on 
crystallization was evaluated by performing both sim- 



ulations including and excluding His, by employing a 
fluctuating lattice Boltzmann method 6] and Langevin 
dynamics [12J, respectively. We prepare our systems 
as an undercooled liquid and let the system crystallize. 
Due to the presence of the confining walls, the nucle- 
ation barrier is sufficiently small, so that we do not 
need special rare event sampling techniques. All simu- 
lations were performed using the MD simulation package 
ESPResSo BHIH. 

As interparticle potential we used a screened Coulomb 
interaction potential 



U(r) = l B k B T 
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where l B is the Bjerrum length, k B is the Boltzmann 
constant, Q\ and Qi give the charges on the interact- 
ing particles, r is distance between the particles. The 
range of the potential is determined by the Debye-Huckcl 
screening length Aj> The static properties of a Yukawa 
system can be characterized by two dimensionless param- 
eters [lfj 
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where w = (3/(4np)) 1 ' 3 is the Wigner-Seitz radius of 
the crystal phase and p the particle density. Phase di- 
agrams of systems with Yukawa-type interactions have 



been calculated both b 
and MD simulations 



by Moi 

Go, EH, 



Monte Carlo simulations [21 1 
which consistently found 
three regimes: a fluid phase and two different solid phases 
with BCC or FCC structure, respectively. For our simu- 
lations, we chose K = 3.0 and V =1260, which is slightly 
above the fluid-solid transition line in the BCC regime. 
The walls act on the particles via a Weeks-Chandler- 
Andersen (WCA) potential. When modeling the solvent 
by a Langevin thermostat [3(, drag and random forces 



on the particles lead to the correct thermal distribution, 
but hydrodynamic interactions are suppressed. The only 
tunable parameter is the friction 7, which is inversely 
proportional to the single particle diffusion constant. 

In order to introduce hydrodynamic interactions with- 
out simulating the solvent explicitly, we used the lattice 
Boltzmann method [6] on a three-dimensional (3D) lat- 
tice with 19 velocity densities (D3Q19). In all reported 
simulations with His, we used a grid spacing of g = 1.0. 
The NVT ensemble was realized by the fluctuating LB 
algorithm [l|, |29j . We treated the colloidal particles as 
point particles that are coupled to the LB fluid via a 
friction term with an adjustable friction constant 7 [2j. 
Note that this friction constant is equivalent to the fric- 
tion in the Langevin dynamics, which is why we use the 
same label 7. The no-slip boundaries modeling the fluid- 
wall interaction were realized by the link bounce back 
rule [33j . The strength of the hydrodynamic interactions 
can be quantified by the hydrodynamic radius 
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where vh depends on the viscosity 77 of the fluid as well 
as of the single particle diffusion coefficient D°, which 
is inversely proportional to the friction constant 7. In 
contrast to other methods for including His, such as 
DPD [8J, |l3[ , LB methods allow the friction 7 to be tuned 
independently from the viscosity rj of the fluid. The par- 
ticle mobility is not simply the inverse of 7, but also 
depends on the viscosity r\ and the lattice spacingof the 
LB grid due to feedback from the moving fluid @, [31| ■ 
Also, the back flow of the solvent introduces finite-size 
effects in a system with periodic boundary conditions. 

The equations of motion of the Yukawa particles were 
integrated by a Velocity- Verlet integrator 12 1. If not 



otherwise stated, the simulations were performed with 
16, 384 particles in a box of size 66 x 16 x 16 confined 
by two planar walls located at x = 0.5 and x = 65.5. 
The time step of the Velocity Verlet integrator was set 
to dt = 0.01, the total simulation length 750,000 time 
steps. The same time step was also used for the LB fluid 
update, when applied. As basic length we use the mean 
particle distance in the crystal phase a — 1.1. The time 
step is dr* = O.Olr, with r = ayk^T '/m p , where m p is 
the mass of the colloids. The viscosity is r\ = 0.8 and the 
density of the fluid is pji = 1.0 as well as the particle 
density p = 1.0. The friction 7 varies between 0.5 and 
12.5. 

The effective diffusion constant depends not only on 
the applied thermostat, but also on the interactions be- 
tween the particles. In order to set up comparable sim- 
ulations, we therefore matched the tracer diffusion coef- 
ficient in the bulk liquid. This was done by measuring 
the MSD of tracer particles in pure bulk systems, both 
with and without His. These measurements were done 
in 3D periodic systems consisting of 16, 384 particles in a 
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FIG. 1: The diffusion coefficient of the tracer particles in the 
bulk Df as a function of the single particle diffusion coeffi- 
cients D° . (in simulation units) . Red squares show results for 
the system without His and the blue triangles for the systems 
with HI. The gray dashed lines are a guide to eye. The black 
lines illustrate the matching of the diffusion coefficient. Two 
different single particle diffusion D° in the Langevin dynam- 
ics and the LB coupling result in the same long-time diffusion 
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box of size 64 x 16 x 16 without confining walls. Figure Q] 
shows the measured long-time diffusion coefficient Df of 
tracer particles in the bulk calculated from the MSD, as 
a function of the single particle diffusion coefficient D°. 
The simulations with and without His led to significantly 
different tracer diffusion coefficients in the bulk, due to 
the fact that e.g. the single particle diffusion in case of 
LB is not simply D° oc 1/7 [2:]. Furthermore, there is 
no simple analytical rule to calculate Df in case of a 
soft-sphere bulk systems. In order to match the tracer 
particle diffusion coefficient Df, we used different val- 
ues of 7 in the Langevin dynamics and the LB coupling. 
The black line in Fig. [T] illustrates this matching: in or- 
der to obtain the same diffusion coefficient Df = 0.015, 
we have to apply 7 = 4.0 for the Langevin dynamics and 
7 = 7.0 for the LB coupling. In the following, we always 
report data with matched tracer diffusion coefficients in 
the bulk, where the given values of 7 are always the ones 
that apply to the LB coupling. 

Using the matched tracer diffusion, we investigated 
the freezing of the undercooled fluid confined between 
two planar walls. In order to distinguish the liquid and 
the different solid phases, we used the Steinhardt or- 
der parameter [30j. Investigations by Moroni et al. [22j 
showed that especially 54 and q 6 are good choices to 
determine whether cubic or hexagonal structures are 
present in the system, respectively. In the following, 
we will focus on FCC and BCC crystal structures, for 
which the qe order parameter is well suited. Due to 



the strong fluctuations in our system, we applied an en- 
hanced averaging method for the Steinhardt order pa- 
rameter q e , introduced by Lecher [16|. The literature 
values are g 6 (BCC) = 0.408018, g 6 (HCP) = 0.42181 and 
g 6 (LIQ) = 0.161962. 

Figure [5] shows the measured q 6 of three snapshots 
taken at different times during a typical simulation run. 
Note that the points only represent the peaks the distri- 
bution of q 6 , since in the crystal, there is a strong layering 
parallel to the wall. In between the peaks, the density 
drops nearly to zero in the crystal, and consequently so 
does the order parameter. As expected, the crystal starts 
with a HCP wall layer, followed by a BCC crystal front 
that grows with time. To evaluate the position of the 
crystal front, we fitted the q 6 peaks to a function of shape 
—h ■ arctan((a; — s)/w), where x is the x-position in the 
simulation box, h is the height difference between q 6 in 
the liquid and in the FCC phase, w is the width of the 
liquid-crystal transition region and s is the position of 
the crystal front. Note that q 6 in the liquid bulk is larger 
than the literature value, since we report only the peaks, 
not the usual average value. 
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FIG. 2: The figure shows our fitting procedure for s (position 
of the crystal front), from the function of q 6 depending on 
x* = x/a for different times fa < fa < fa. The gray dashed 
lines indicate the computed front locations (si < S2 < S3), 
while the black dashed line shows the fit for si at time fa. 
The symbols represent the peaks of the q 6 order parameter. 



After some initial time, the crystal grew very uni- 
formly, so that we could determine a constant growth 
velocity u by a linear fitting of the front position. Fig- 
ure [3] shows the measured velocities u as a function of the 
hydrodynamic radius r#, which we varied by changing 
the friction coefficient 7 and applying the matching pro- 
cedure described above. Every measurement represents 
the mean growth velocity sampled from 24 independent 
runs. For hydrodynamic radii r* H = ruj a < 0.025, where 
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FIG. 3: Growth velocity u normalized by the bulk diffusion 
coefficient D t times the mean particle distance a as a function 
of the relative hydrodynamic radius r* H = th/cl. The red 
squares show the results for simulations without His, which 
are almost independent of r* H within the error bars. The blue 
triangles represent the results for simulations with His, which 
show a strong decay of the growth velocity with increasing 
hydrodynamic radius. The gray dashed lines are guidelines 
to eye. 



a is the mean particle distance in the crystal phase, the 
influence of HI is almost negligible as one would expect. 
But already in case of moderate ratios 0.1 < r* H < 0.25, 
hydrodynamic interactions reduce the crystal growth ve- 
locity by up to a factor of 3 at r* H — 0.25. In case of no 
His the normalized growth velocity is virtually constant, 
with a decay for small frictions due to improper coupling 
to the thermostat. 

In order to elucidate what causes this difference, we 
analyzed the particle diffusion in the system relative to 
the actual position of the crystal front. To accomplish 
this, we binned our system along the growth direction 
into bins of width b = 2a and determined the long-time 
diffusion coefficient D^ OM of the center of mass in the 
direction of growth, which can be seen as a measure for 
the transport of particles towards the growing crystal 
front. In Fig. Q] D^ OM relative to the position of the 
crystal front x rel , normalized by the center of mass dif- 
fusion in the bulk D2P of the Langevin simulation is 



shown. The front of the crystal is located at x reL = 0, 
while the pure bulk fluid phase is located at x rel — 7 and 
the crystal phase is present for x rel < 0. As expected, 
the long-time diffusion coefficients for the center of mass 
in the crystalline region are almost zero, rise in the re- 
gion of the crystal front, and settle off to the liquid bulk 
value far away of the crystal front. The left-hand side of 
Fig. Q] shows the values for low r* H = rn/a = 0.025 ratio, 
which are virtually the same for both systems: with and 
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FIG. 4: Long-time diffusion coefficient of the center of mass 
in the direction of growth. Every value is the long-time diffu- 
sion coefficient of the center of mass for particles in a bin of 
size b — 2a normalized by the value of the Langevin simula- 
tion (w/o His). The front of the crystal is located at x rel = 0, 
while x rel — 7 is located in the pure bulk fluid phase, (a) The 
first two cases (I: red circles, II: blue rhombus) have a low 
r* H = th /a = 0.025 ratio and show virtually similar behavior, 
(b) For moderate r* H = 0.25, case IV (blue triangles) with His 
shows a different behavior from case III (red squares) without 
His, especially in the region in front of the crystal phase. 



without His. However, on the right-hand side r* H = 0.25 
is shown, corresponding to moderate hydrodynamic cou- 
pling, for which a significantly reduction of the diffusion 
is present in case of His. Especially, nearby the crystal 
front a drop down of more than a factor of two in the dif- 
fusion coefficient occurs, corresponding to a drastically 
reduced transport towards the crystal. We believe that 
the reduced diffusion towards the crystal front in the LB 
system is caused by hydrodynamic screening due to the 
presence of the wall-like crystal surface. Hydrodynamic 
screening nearby walls and their influence on the diffu- 
sion coefficients parallel and perpendicular to the wall 
has been reported by Feitosa et al. [9j and von Hansen et 



al. [32 1 . 



Our simulations show that hydrodynamic interactions 
have a strong influence on crystallization, even at mod- 
erate hydrodynamic radii. Similar finding has been re- 
ported by Schilling et al. 27] as well. The effects arise 



mainly on the particle transport towards the crystal 
front, which are in particular important for nucleation 
processes. This puts the common assumption into doubt 
that hydrodynamic interactions can be ignored when 
studying crystallization or nucleation in suspensions. At 
least in Yukawa suspensions, these processes do not seem 
to be quasi-static, and the often drawn analogy to true 
melts might not be true. 
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